/*
 * case.cpp
 *
 *  Created on: 11-09-2012
 *      Author: mith
 */

#include "case_2D2t_SGA_fast.h"
// MY HEADERS
#include"mros-v42.h"

namespace case_2D2t_SGA_fast
{
global glob;

// *******************************************************************************************  //
// *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY ***
// *******************************************************************************************  //

#define Sqrt sqrt
#define Power pow

/*
 *********************************************************
 ******************** Structures *************************
 *********************************************************
*/

struct parameters {
	double n;
	double J;
	double t;
	double tPrim;
	double U;
	double beta;

	unsigned int L; //8
};

struct variables {
	double barchiAB;
	double chiS;
	double chiT;
	double DeltaS;
	double DeltaT;
	double mu;
	double mAF;
	double d;
    double Energy; //13
 };

/*
 *********************************************************
 ***************** Small useful functions ****************
 *********************************************************
*/

inline void clear(parameters &p) {
	p.n = p.J = p.t = p.tPrim = p.U = p.beta = 0.; p.L = 0;
}

inline void clear(variables &v) {
	v.barchiAB = v.chiS = v.chiT = v.DeltaS = v.DeltaT = v.mu = v.mAF = v.d = 0;
	v.Energy = 0;
}

inline void addTo(const variables &q, variables &v, double d=1.) {
	v.barchiAB += q.barchiAB*d;
	v.chiS += q.chiS*d;
	v.chiT += q.chiT*d;
	v.DeltaS += q.DeltaS*d;
	v.DeltaT += q.DeltaT*d;
	v.mu += q.mu*d;
	v.mAF += q.mAF*d;
	v.d += q.d*d;
	v.Energy += q.Energy*d;
}

inline void copyTo(variables &v, const variables &q, double d=1.) {
	v.barchiAB = q.barchiAB*d;
	v.chiS = q.chiS*d;
	v.chiT = q.chiT*d;
	v.DeltaS = q.DeltaS*d;
	v.DeltaT = q.DeltaT*d;
	v.mu = q.mu*d;
	v.mAF = q.mAF*d;
	v.d = q.d*d;
	v.Energy = q.Energy*d;
}

inline void copyTo(parameters &v, const parameters &q) {
	v.n = q.n;
	v.J = q.J;
	v.t = q.t;
	v.tPrim = q.tPrim;
	v.U = q.U;
	v.beta = q.beta;
	v.L = q.L;
}

inline double abs(double a) {
	if(a<0) a=-a;
	return a;
}

/*
 *********************************************************
 ***************** CLASS Functions ***********************
 *********************************************************
*/

class Func {
private:
	parameters p; // parameters  p.n  p.J  p.t  p.tPrim  p.U  p.beta  p.L
	variables v; // variables  v.chiAB  v.chiS  v.chiT  v.barDelta  v.simDelta  v.mu  v.mAF  v.d
	double n, d, m, mu; // contractions

	// Abbreviations - dependend on kx, ky
	double gammak;
	double etak;
	double gammaPrimk;

public:
	double E[4];
	double dEdlbarchiAB[4];
	double dEdlchiS[4];
	double dEdlchiT[4];
	double dEdlDeltaS[4];
	double dEdlDeltaT[4];
	double dEdln[4];
	double dEdlmAF[4];

	double gt, gs, gtPrimMin, gtPrimMax, gtPrim2Min, gtPrim2Max;
	 double dgtdm, dgtdd, dgtdn;
	 double dgsdm, dgsdd, dgsdn;
	 double dgtPrimMindm, dgtPrimMindd, dgtPrimMindn;
	 double dgtPrimMaxdm, dgtPrimMaxdd, dgtPrimMaxdn;
	double lbarchiAB, lsimchiAB, lchiS, lchiT, lchiSPrim, lchiTPrim, lDeltaS, lDeltaT, lsimDeltaB, ln, lmAF;

	Func(const parameters &p, const variables &v) {
		set(p, v);
		//refresh();
	}

	void set(const parameters &p, const variables &v) {
		copyTo(this->p, p);
		copyTo(this->v, v);
		refresh();
	}

	void set(const parameters &p) {
		copyTo(this->p, p);
		refresh();
	}

	// Gutzwiller
	double fgt();
	double fgs();
	double fgBDelta();
	double fgm();
	double fgtPrimMin();
	double fgtPrimMax();

	double fdgtdm();
	double fdgtdd();
	double fdgtdn();

	double fdgsdm();
	double fdgsdd();
	double fdgsdn();

	double fdgtPrimMindm();
	double fdgtPrimMindd();
	double fdgtPrimMindn();
	double fdgtPrimMaxdm();
	double fdgtPrimMaxdd();
	double fdgtPrimMaxdn();

    // Energy
    double E1(double gammak, double gammaPrimk, double etak);
    double E3(double gammak, double gammaPrimk, double etak);

    double dE1dl_barchiAB(double gammak, double gammaPrimk, double etak);
    double dE1dl_chiS(double gammak, double gammaPrimk, double etak);
    double dE1dl_chiT(double gammak, double gammaPrimk, double etak);
    double dE1dl_n(double gammak, double gammaPrimk, double etak);
    double dE1dl_mAF(double gammak, double gammaPrimk, double etak);
    double dE1dl_DeltaS(double gammak, double gammaPrimk, double etak);
    double dE1dl_DeltaT(double gammak, double gammaPrimk, double etak);

    double dE3dl_barchiAB(double gammak, double gammaPrimk, double etak);
    double dE3dl_chiS(double gammak, double gammaPrimk, double etak);
    double dE3dl_chiT(double gammak, double gammaPrimk, double etak);
    double dE3dl_n(double gammak, double gammaPrimk, double etak);
    double dE3dl_mAF(double gammak, double gammaPrimk, double etak);
    double dE3dl_DeltaS(double gammak, double gammaPrimk, double etak);
    double dE3dl_DeltaT(double gammak, double gammaPrimk, double etak);

	// DEBUG
	void printLambda() {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e (l)\n", lbarchiAB, lchiS, lchiT, lDeltaS, lDeltaT, ln, lmAF);
		printf("### %25.15e %25.15e %25.15e %25.15e (g)\n", gt, gs, gtPrimMin, gtPrimMax);
	}

	// REFRESH
	void refresh() {
		n = p.n;
		d = v.d;
		mu = v.mu;
		m = v.mAF;

		// obliczam funkcje Gutzwillara i poszczegolne pochodne
		gt = fgt();
		gs = fgs();
		gtPrimMin = fgtPrimMin();
		gtPrimMax = fgtPrimMax();

		dgtdn = fdgtdn();
		dgtdd = fdgtdd();
		dgtdm = fdgtdm();

		dgsdn = fdgsdn();
		dgsdm = fdgsdm();
		dgsdd = fdgsdd();

		dgtPrimMindm = fdgtPrimMindm();
		dgtPrimMindd = fdgtPrimMindd();
		dgtPrimMindn = fdgtPrimMindn();
		dgtPrimMaxdm = fdgtPrimMaxdm();
		dgtPrimMaxdd = fdgtPrimMaxdd();
		dgtPrimMaxdn = fdgtPrimMaxdn();


		// obliczanie poszczegolnych lambda:
		lbarchiAB = -gt*p.t + (3./4.)*gs*p.J*v.barchiAB;
		lchiS = -(1./2.)*p.tPrim*(gtPrimMax + gtPrimMin);
		lchiT = 0.;
		lDeltaS = (3./4.)*gs*p.J*v.DeltaS;
		lDeltaT = -(1./4.)*gs*p.J*v.DeltaT;
		ln = -16.*p.t*v.barchiAB*dgtdn \
			- 8.*p.tPrim*v.chiS*dgtPrimMindn - 8.*p.tPrim*v.chiS*dgtPrimMaxdn \
			- 2.*p.J*(-0.5*v.mAF*v.mAF - 3.*v.barchiAB*v.barchiAB -3.*v.DeltaS*v.DeltaS + v.DeltaT*v.DeltaT)*dgsdn;
		lmAF = 2.*gs*p.J*v.mAF - 16.*p.t*v.barchiAB*dgtdm \
			- 8.*p.tPrim*v.chiS*dgtPrimMindm - 8.*p.tPrim*v.chiS*dgtPrimMaxdm \
			- 2.*p.J*(-0.5*v.mAF*v.mAF - 3.*v.barchiAB*v.barchiAB -3.*v.DeltaS*v.DeltaS + v.DeltaT*v.DeltaT)*dgsdm;

		// zeby pozniej mozna bylo wypisac
		glob.lbarchiAB = lbarchiAB;
		glob.lchiS = lchiS;
		glob.lchiT = lchiT;
		glob.lDeltaS = lDeltaS;
		glob.lDeltaT = lDeltaT;
		glob.ln = ln;
		glob.lmAF = lmAF;
	}

	//REFRESH XY
	void refresh(double kx, double ky) {
		gammak = 2.*(cos(kx) + cos(ky));
		gammaPrimk = 4.*cos(kx)*cos(ky);
		etak = 2.*(cos(kx) - cos(ky));

		refresh();

		E[0]          		= E1(gammak, gammaPrimk, etak);
		dEdlbarchiAB[0]  	= dE1dl_barchiAB(gammak, gammaPrimk, etak);
		dEdlchiS[0] 		= dE1dl_chiS(gammak, gammaPrimk, etak);
		dEdlchiT[0] 		= dE1dl_chiT(gammak, gammaPrimk, etak);
		dEdlDeltaS[0] 		= dE1dl_DeltaS(gammak, gammaPrimk, etak);
		dEdlDeltaT[0] 		= dE1dl_DeltaT(gammak, gammaPrimk, etak);
		dEdln[0]      		= dE1dl_n(gammak, gammaPrimk, etak);
		dEdlmAF[0]    		= dE1dl_mAF(gammak, gammaPrimk, etak);

		E[1]          		= -E[0];
		dEdlbarchiAB[1]  	= -dEdlbarchiAB[0];
		dEdlchiS[1] 		= -dEdlchiS[0];
		dEdlchiT[1] 		= -dEdlchiT[0];
		dEdlDeltaS[1] 		= -dEdlDeltaS[0];
		dEdlDeltaT[1] 		= -dEdlDeltaT[0];
		dEdln[1]      		= -dEdln[0];
		dEdlmAF[1]    		= -dEdlmAF[0];

		E[2]          		= E3(gammak, gammaPrimk, etak);
		dEdlbarchiAB[2]  	= dE3dl_barchiAB(gammak, gammaPrimk, etak);
		dEdlchiS[2] 		= dE3dl_chiS(gammak, gammaPrimk, etak);
		dEdlchiT[2] 		= dE3dl_chiT(gammak, gammaPrimk, etak);
		dEdlDeltaS[2] 		= dE3dl_DeltaS(gammak, gammaPrimk, etak);
		dEdlDeltaT[2] 		= dE3dl_DeltaT(gammak, gammaPrimk, etak);
		dEdln[2]      		= dE3dl_n(gammak, gammaPrimk, etak);
		dEdlmAF[2]    		= dE3dl_mAF(gammak, gammaPrimk, etak);

		E[3]          		= -E[2];
		dEdlbarchiAB[3]  	= -dEdlbarchiAB[2];
		dEdlchiS[3] 		= -dEdlchiS[2];
		dEdlchiT[3] 		= -dEdlchiT[2];
		dEdlDeltaS[3] 		= -dEdlDeltaS[2];
		dEdlDeltaT[3] 		= -dEdlDeltaT[2];
		dEdln[3]      		= -dEdln[2];
		dEdlmAF[3]    		= -dEdlmAF[2];
	}
};


/////////////////
///////////////// FullSimplify
/////////////////

double Func::fgBDelta() {
    //n = 1.-p.sdelta;
    double r = n/2. + v.mAF/2.;
    double w = n/2. - v.mAF/2.;
    double d=v.d*v.d;
	double part1 = (n-2.*d)/(2.*(n-2.*r*w));
	double part2 = sqrt((1.-w)*(1.-n+d)/(1.-r)) + sqrt(w*d/r);
	double part3 = sqrt((1.-r)*(1.-n+d)/(1.-w)) + sqrt(r*d/w);
	return part1*(part2*part2 + part3*part3);
}

double Func::fgm() {
   double r = n/2. + v.mAF/2.;
   double w = n/2. - v.mAF/2.;
   double d = v.d*v.d;
   return (n-2.*d)/(n-2.*w*r);
}

double Func::fgt() {
	return (2*(2*pow(d,2) - n)*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n))))/
			   (-pow(m,2) + (-2 + n)*n);
}

double Func::fgs() {
	return (4*pow(-2*pow(d,2) + n,2))/pow(pow(m,2) - (-2 + n)*n,2);
}

double Func::fgtPrimMin() {
	return (2*(2*pow(d,2) - n)*pow(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n)),2))/(-pow(m,2) + (-2 + n)*n);
}

double Func::fgtPrimMax() {
	return (2*(2*pow(d,2) - n)*pow(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)),2))/(-pow(m,2) + (-2 + n)*n);
}

////////////////////////////////

double Func::fdgtdm() {
	return (4*(-2*pow(d,2) + n)*(-(m*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*
	          (sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n)))) +
	       (pow(m,4)*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n))*sqrt((pow(d,2)*(-m + n))/(m + n)) - sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n))*sqrt((pow(d,2)*(m + n))/(-m + n))))/
	        ((m - n)*(2 + m - n)*(-2 + m + n)*(m + n)) - (pow(-2 + n,2)*pow(n,2)*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n))*sqrt((pow(d,2)*(-m + n))/(m + n)) -
	            sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n))*sqrt((pow(d,2)*(m + n))/(-m + n))))/((m - n)*(2 + m - n)*(-2 + m + n)*(m + n))))/pow(pow(m,2) - (-2 + n)*n,2);
}

double Func::fdgtdd() {
	return (2*(-4*d*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n))) +
		       (-2*pow(d,2) + n)*((d*sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)))/(1 + pow(d,2) - n) + sqrt((pow(d,2)*(-m + n))/(m + n))/d)*
		        (sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n))) +
		       (-2*pow(d,2) + n)*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*
		        ((d*sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)))/(1 + pow(d,2) - n) + sqrt((pow(d,2)*(m + n))/(-m + n))/d)))/(pow(m,2) - (-2 + n)*n);
}

double Func::fdgtdn() {
	return (4*(((pow(m,2) - (-2 + n)*n)*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) +
            sqrt((pow(d,2)*(-m + n))/(m + n)))*
          (sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt(-((pow(d,2)*(m + n))/(m - n)))))/2.
         - (1 - n)*(-2*pow(d,2) + n)*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) +
          sqrt((pow(d,2)*(-m + n))/(m + n)))*
        (sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n))) +
       ((-2*pow(d,2) + n)*(pow(m,2) - (-2 + n)*n)*
          ((pow(m,2) + 2*m*(1 + pow(d,2) - n) - pow(-2 + n,2))/
             (sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n))*pow(-2 + m + n,2)) +
            (2*pow(d,2)*m)/(sqrt((pow(d,2)*(-m + n))/(m + n))*pow(m + n,2)))*
          (sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n))))/4.\
        + ((-2*pow(d,2) + n)*(pow(m,2) - (-2 + n)*n)*
          (sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*
          ((pow(m,2) - 2*m*(1 + pow(d,2) - n) - pow(-2 + n,2))/
             (pow(2 + m - n,2)*sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n))) +
            (2*m*sqrt(-((pow(d,2)*(m + n))/(m - n))))/(pow(m,2) - pow(n,2))))/4.))/
   pow(pow(m,2) - (-2 + n)*n,2);
}

double Func::fdgsdm() {
	return (-16*m*pow(-2*pow(d,2) + n,2))/pow(pow(m,2) - (-2 + n)*n,3);
}

double Func::fdgsdd() {
	return (64*pow(d,3) - 32*d*n)/pow(pow(m,2) - (-2 + n)*n,2);
}

double Func::fdgsdn() {
	return (8*(-2*pow(d,2) + n)*(pow(m,2) - 4*pow(d,2)*(-1 + n) + pow(n,2)))/pow(pow(m,2) - (-2 + n)*n,3);
}

double Func::fdgtPrimMindm() {
	return (4*(-2*pow(d,2) + n)*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n)))*
		     (-(m*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n)))) +
		       ((pow(m,2) - (-2 + n)*n)*((2*(-2 + n)*sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)))/((-2 - m + n)*(-2 + m + n)) + (2*n*sqrt((pow(d,2)*(m + n))/(-m + n)))/(-pow(m,2) + pow(n,2))))/2.))/
		   pow(pow(m,2) - (-2 + n)*n,2);
}

double Func::fdgtPrimMindd() {
	return (4*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n)))*
		     (-2*d*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n))) +
		       (-2*pow(d,2) + n)*((d*sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)))/(1 + pow(d,2) - n) + sqrt((pow(d,2)*(m + n))/(-m + n))/d)))/(pow(m,2) - (-2 + n)*n);
}

double Func::fdgtPrimMindn() {
	return (4*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n)))*
		     (((pow(m,2) - (-2 + n)*n)*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt(-((pow(d,2)*(m + n))/(m - n)))))/2. +
		       (-1 + n)*(-2*pow(d,2) + n)*(sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n)) + sqrt((pow(d,2)*(m + n))/(-m + n))) +
		       ((-2*pow(d,2) + n)*(pow(m,2) - (-2 + n)*n)*((pow(m,2) - 2*m*(1 + pow(d,2) - n) - pow(-2 + n,2))/(pow(2 + m - n,2)*sqrt(((-1 - pow(d,2) + n)*(-2 + m + n))/(2 + m - n))) +
		            (2*m*sqrt(-((pow(d,2)*(m + n))/(m - n))))/(pow(m,2) - pow(n,2))))/2.))/pow(pow(m,2) - (-2 + n)*n,2);
}

double Func::fdgtPrimMaxdm() {
	return (4*(-2*pow(d,2) + n)*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*
		     (-(m*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))) +
		       ((pow(m,2) - (-2 + n)*n)*((2*(-2 + n)*sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)))/((2 + m - n)*(-2 + m + n)) - (2*pow(d,2)*n)/(sqrt((pow(d,2)*(-m + n))/(m + n))*pow(m + n,2))))/2.))/
		   pow(pow(m,2) - (-2 + n)*n,2);
}

double Func::fdgtPrimMaxdd() {
	return (4*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*
		     (-2*d*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n))) +
		       (-2*pow(d,2) + n)*((d*sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)))/(1 + pow(d,2) - n) + sqrt((pow(d,2)*(-m + n))/(m + n))/d)))/(pow(m,2) - (-2 + n)*n);
}

double Func::fdgtPrimMaxdn() {
	return (4*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n)))*
		     ((-1 + n)*(-2*pow(d,2) + n)*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n))) +
		       ((pow(m,2) - (-2 + n)*n)*(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n)) + sqrt((pow(d,2)*(-m + n))/(m + n))))/2. +
		       ((-2*pow(d,2) + n)*(pow(m,2) - (-2 + n)*n)*((pow(m,2) + 2*m*(1 + pow(d,2) - n) - pow(-2 + n,2))/(sqrt(((2 + m - n)*(-1 - pow(d,2) + n))/(-2 + m + n))*pow(-2 + m + n,2)) +
		            (2*pow(d,2)*m)/(sqrt((pow(d,2)*(-m + n))/(m + n))*pow(m + n,2))))/2.))/pow(pow(m,2) - (-2 + n)*n,2);
}

/////////////////////////////////////////////////////////

double Func::E1(double gammak, double gammaPrimk, double etak) {
	return -Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
		      Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) -
		      2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		         4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2))))/2.;
}

double Func::E3(double gammak, double gammaPrimk, double etak) {
	return -Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
		      Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) +
		      2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		         4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2))))/2.;
}

double Func::dE1dl_barchiAB(double gammak, double gammaPrimk, double etak) {
	return (2*Power(gammak,2)*lbarchiAB*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2) -
		       Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		         4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))))/
		   (Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		       4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))*
		     Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
		       Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) -
		       2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))));
}

double Func::dE3dl_barchiAB(double gammak, double gammaPrimk, double etak) {
	return (-2*Power(gammak,2)*lbarchiAB*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2) +
		       Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		         4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))))/
		   (Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		       4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))*
		     Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
		       Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) +
		       2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))));
}

double Func::dE1dl_chiS(double gammak, double gammaPrimk, double etak) {
	return -((gammaPrimk*(2*gammaPrimk*lchiS + ln + 2*mu + (-4*Power(gammak,2)*Power(lbarchiAB,2)*(2*gammaPrimk*lchiS + ln + 2*mu) -
            (2*gammaPrimk*lchiT + lmAF)*(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu)))/
          Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
            4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))))/
     Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
       Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) -
       2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))));
}

double Func::dE3dl_chiS(double gammak, double gammaPrimk, double etak) {
	return -((gammaPrimk*(2*gammaPrimk*lchiS + ln + 2*mu + (4*Power(gammak,2)*Power(lbarchiAB,2)*(2*gammaPrimk*lchiS + ln + 2*mu) +
            (2*gammaPrimk*lchiT + lmAF)*(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu)))/
          Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
            4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))))/
     Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
       Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) +
       2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))));
}

double Func::dE1dl_chiT(double gammak, double gammaPrimk, double etak) {
	return (gammaPrimk*(8*Power(gammaPrimk,3)*Power(lchiS,2)*lchiT + 4*Power(etak,2)*lDeltaS*lDeltaT*(ln + 2*mu) + lmAF*Power(ln + 2*mu,2) +
		       4*Power(gammaPrimk,2)*lchiS*(lchiS*lmAF + 2*lchiT*(ln + 2*mu)) -
		       lmAF*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2))) +
		       2*gammaPrimk*(4*Power(etak,2)*lchiS*lDeltaS*lDeltaT + (ln + 2*mu)*(2*lchiS*lmAF + lchiT*(ln + 2*mu)) -
		          lchiT*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		             4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2))))))/
		   (Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		       4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))*
		     Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
		       Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) -
		       2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))));
}

double Func::dE3dl_chiT(double gammak, double gammaPrimk, double etak) {
	return (gammaPrimk*(-8*Power(gammaPrimk,3)*Power(lchiS,2)*lchiT - 4*Power(etak,2)*lDeltaS*lDeltaT*(ln + 2*mu) - lmAF*Power(ln + 2*mu,2) -
		       4*Power(gammaPrimk,2)*lchiS*(lchiS*lmAF + 2*lchiT*(ln + 2*mu)) -
		       lmAF*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2))) -
		       2*gammaPrimk*(4*Power(etak,2)*lchiS*lDeltaS*lDeltaT + (ln + 2*mu)*(2*lchiS*lmAF + lchiT*(ln + 2*mu)) +
		          lchiT*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		             4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2))))))/
		   (Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		       4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))*
		     Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
		       Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) +
		       2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))));
}

double Func::dE1dl_DeltaS(double gammak, double gammaPrimk, double etak) {
	return (2*Power(etak,2)*(4*Power(etak,2)*lDeltaS*Power(lDeltaT,2) + lDeltaT*(2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu) -
		       lDeltaS*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))))/
		   (Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		       4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))*
		     Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
		       Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) -
		       2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))));
}

double Func::dE3dl_DeltaS(double gammak, double gammaPrimk, double etak) {
	return (2*Power(etak,2)*(-4*Power(etak,2)*lDeltaS*Power(lDeltaT,2) - lDeltaT*(2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu) -
		       lDeltaS*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))))/
		   (Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		       4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))*
		     Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
		       Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) +
		       2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))));
}

double Func::dE1dl_DeltaT(double gammak, double gammaPrimk, double etak) {
	return (2*Power(etak,2)*(4*Power(gammaPrimk,2)*lchiS*lchiT*lDeltaS + 4*Power(gammak,2)*Power(lbarchiAB,2)*lDeltaT + 2*gammaPrimk*lDeltaS*(lchiS*lmAF + lchiT*(ln + 2*mu)) +
		       lDeltaS*(4*Power(etak,2)*lDeltaS*lDeltaT + lmAF*(ln + 2*mu)) - lDeltaT*
		        Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))))/
		   (Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		       4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))*
		     Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
		       Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) -
		       2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))));
}

double Func::dE3dl_DeltaT(double gammak, double gammaPrimk, double etak) {
	return (-2*Power(etak,2)*(4*Power(gammaPrimk,2)*lchiS*lchiT*lDeltaS + 4*Power(gammak,2)*Power(lbarchiAB,2)*lDeltaT + 2*gammaPrimk*lDeltaS*(lchiS*lmAF + lchiT*(ln + 2*mu)) +
		       lDeltaS*(4*Power(etak,2)*lDeltaS*lDeltaT + lmAF*(ln + 2*mu)) + lDeltaT*
		        Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))))/
		   (Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		       4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))*
		     Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
		       Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) +
		       2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))));
}


double Func::dE1dl_n(double gammak, double gammaPrimk, double etak) {
	return -((gammaPrimk*lchiS + ln/2. + mu + (-4*Power(gammak,2)*Power(lbarchiAB,2)*(2*gammaPrimk*lchiS + ln + 2*mu) -
	          (2*gammaPrimk*lchiT + lmAF)*(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu)))/
	        (2.*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
	            4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))))/
	     Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
	       Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) -
	       2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
	          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))));
}

double Func::dE3dl_n(double gammak, double gammaPrimk, double etak) {
	return -((gammaPrimk*lchiS + ln/2. + mu + (4*Power(gammak,2)*Power(lbarchiAB,2)*(2*gammaPrimk*lchiS + ln + 2*mu) +
	          (2*gammaPrimk*lchiT + lmAF)*(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu)))/
	        (2.*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
	            4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))))/
	     Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
	       Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) +
	       2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
	          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))));
}

double Func::dE1dl_mAF(double gammak, double gammaPrimk, double etak) {
	return (8*Power(gammaPrimk,3)*Power(lchiS,2)*lchiT + 4*Power(etak,2)*lDeltaS*lDeltaT*(ln + 2*mu) + lmAF*Power(ln + 2*mu,2) +
		     4*Power(gammaPrimk,2)*lchiS*(lchiS*lmAF + 2*lchiT*(ln + 2*mu)) - lmAF*
		      Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		        4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2))) +
		     2*gammaPrimk*(4*Power(etak,2)*lchiS*lDeltaS*lDeltaT + (ln + 2*mu)*(2*lchiS*lmAF + lchiT*(ln + 2*mu)) -
		        lchiT*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		           4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))))/
		   (2.*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		       4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))*
		     Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
		       Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) -
		       2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))));
}

double Func::dE3dl_mAF(double gammak, double gammaPrimk, double etak) {
	return (-8*Power(gammaPrimk,3)*Power(lchiS,2)*lchiT - 4*Power(etak,2)*lDeltaS*lDeltaT*(ln + 2*mu) - lmAF*Power(ln + 2*mu,2) -
		     4*Power(gammaPrimk,2)*lchiS*(lchiS*lmAF + 2*lchiT*(ln + 2*mu)) - lmAF*
		      Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		        4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2))) -
		     2*gammaPrimk*(4*Power(etak,2)*lchiS*lDeltaS*lDeltaT + (ln + 2*mu)*(2*lchiS*lmAF + lchiT*(ln + 2*mu)) +
		        lchiT*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		           4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))))/
		   (2.*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		       4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))*
		     Sqrt(4*Power(gammak,2)*Power(lbarchiAB,2) + 4*Power(gammaPrimk,2)*(Power(lchiS,2) + Power(lchiT,2)) + 4*Power(etak,2)*(Power(lDeltaS,2) + Power(lDeltaT,2)) + Power(lmAF,2) +
		       Power(ln + 2*mu,2) + 4*gammaPrimk*(lchiT*lmAF + lchiS*(ln + 2*mu)) +
		       2*Sqrt(Power(4*Power(etak,2)*lDeltaS*lDeltaT + (2*gammaPrimk*lchiT + lmAF)*(2*gammaPrimk*lchiS + ln + 2*mu),2) +
		          4*Power(gammak,2)*Power(lbarchiAB,2)*(4*Power(etak,2)*Power(lDeltaT,2) + Power(2*gammaPrimk*lchiS + ln + 2*mu,2)))));
}


/*
 *********************************************************
 ******************* CLASS Equations *********************
 *********************************************************
*/

class Equa {
private:
	parameters p; // parameters
	variables v; // variables
	Func *f;

	void onePoint(variables &y, int i, int j);

public:
	Equa(parameters &p, variables &v) : p(p), v(v) { f = new Func(p, v); }
	~Equa() { delete f; }

	void setPar(parameters &p) {
		this->p = p;
		f->set(p);
	}

	void run(variables &w);
};

void Equa::onePoint(variables &w, int i, int j) {
	double x = (double)i*2.*M_PI/p.L;
	double y = (double)j*2.*M_PI/p.L;

	f->refresh(x, y);
	if(DEEP_DEBUG) { f->printLambda(); }

	double Etan[4] = { (1-tanh(p.beta*f->E[0]/2.))/2., (1-tanh(p.beta*f->E[1]/2.))/2., (1-tanh(p.beta*f->E[2]/2.))/2., (1-tanh(p.beta*f->E[3]/2.))/2. };

	for(int i=0; i<4; i++) {
		w.barchiAB += Etan[i]*f->dEdlbarchiAB[i];
		w.chiS += Etan[i]*f->dEdlchiS[i];
		w.chiT += Etan[i]*f->dEdlchiT[i];
		w.DeltaS += Etan[i]*f->dEdlDeltaS[i];
		w.DeltaT += Etan[i]*f->dEdlDeltaT[i];
		w.mu += Etan[i]*f->dEdln[i];
		w.mAF += Etan[i]*f->dEdlmAF[i];
		// w.d w tym miejscu jeszcze nie

        //if(f->E[i]>0) { w.Energy += -exp(-p.beta*f->E[i])/p.beta; } // basically 0
        //if(f->E[i]<0) { w.Energy += f->E[i]; }

        // it would be safer have following code - warning: 30 is chosen arbitrarily
		if(p.beta*f->E[i] > 30.) { w.Energy += -exp(-p.beta*f->E[i])/p.beta; } // basically 0
		else if(p.beta*f->E[i] < -30) { w.Energy += f->E[i]; }
		else { w.Energy += -log(1 + exp(-p.beta*f->E[i]))/p.beta; }
	}
}

void Equa::run(variables &w)
{
	clear(w);
	variables temp; clear(temp); // temporary vector using during summation
	unsigned int max = p.L/2; // attention: L should be even!

	//int iii = 0;
	// SERIAL CASE
	// SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
	// * edges (2 times)
	for(unsigned int i=1; i<max; i++) {
		onePoint(temp, i, 0);
		onePoint(temp, 0, i);
		onePoint(temp, i, max-i); //iii += 6;
	}
   addTo(temp, w, 2);
   clear(temp);
   // * interior area (4 times) - SERIAL CASE
   for(unsigned int i=1; i<max; i++) {
		for(unsigned int j=1; j<i; j++) {
			onePoint(temp, max-i, j); //iii += 4;
     	}
	}
	addTo(temp, w, 4);
	clear(temp);
	// * missing corners
	onePoint(w, max, 0);
	//onePoint(w, 0, max); // <-- tu bylo o jeden punkt za duzo!
	onePoint(w, 0, 0);
	//iii += 2;
	//printf("iii=%d\n", iii*2);

	// what was out of sumation
	double L2 = (double)(p.L*p.L);
	w.barchiAB /= L2;
	w.chiS /= L2;
	w.chiT /= L2;
	w.DeltaS /= L2;
	w.DeltaT /= L2;
	w.mu /= L2;
	w.mAF /= L2;
    w.Energy /= L2;

	w.barchiAB += 8.*v.barchiAB;
	w.chiS += 8.*v.chiS;
	w.chiT += 8.*v.chiT;
	w.DeltaS += 8.*v.DeltaS;
	w.DeltaT += 8.*v.DeltaT;
	w.mu += -0.5*(1-p.n);
	w.mAF += 0.5*v.mAF;

	// calculate dd
	w.d = 2.*v.d*p.U \
			+ 8.*p.t*v.barchiAB*f->dgtdd \
			+ 4.*p.tPrim*v.chiS*f->dgtPrimMaxdd + 4.*p.tPrim*v.chiS*f->dgtPrimMindd \
			+ p.J*(-0.5*v.mAF*v.mAF - 3.*v.barchiAB*v.barchiAB -3.*v.DeltaS*v.DeltaS + v.DeltaT*v.DeltaT)*f->dgsdd;

    v.Energy = w.Energy \
    		+ 8.*f->gt*p.t*v.barchiAB \
    		+ 4.*f->gtPrimMax*p.tPrim*v.chiS + 4.*f->gtPrimMin*p.tPrim*v.chiS \
    		+ f->gs*p.J*(-0.5*v.mAF*v.mAF - 3.*v.barchiAB*v.barchiAB -3.*v.DeltaS*v.DeltaS + v.DeltaT*v.DeltaT) \
    		+ 8.*(f->lbarchiAB*v.barchiAB + f->lchiS*v.chiS + f->lchiT*v.chiT  + f->lDeltaS*v.DeltaS + f->lDeltaT*v.DeltaT) \
    		+ p.U*v.d*v.d + 0.5*f->ln*(p.n-1) + 0.5*f->lmAF*v.mAF -v.mu;

	//if(w.chiAB != w.chiAB) { w.chiAB = 1000; }
	//if(w.chiS != w.chiS) { w.chiS = 1000; }
	//if(w.chiT != w.chiT) { w.chiT = 1000; }
	//if(w.barDelta != w.barDelta) { w.barDelta = 1000; }
	//if(w.simDelta != w.simDelta) { w.simDelta = 1000; }
	//if(w.mu != w.mu) { w.mu = 1000; }
	//if(w.mAF != w.mAF) { w.mAF = 1000; }

    // results
    case_2D2t_SGA_fast::glob.Energy = v.Energy + v.mu*p.n;
    case_2D2t_SGA_fast::glob.gt = f->fgt();
    case_2D2t_SGA_fast::glob.gs = f->fgs();
    case_2D2t_SGA_fast::glob.gtPrimMax = f->fgtPrimMax();
    case_2D2t_SGA_fast::glob.gtPrimMin = f->fgtPrimMin();
	case_2D2t_SGA_fast::glob.gs = f->fgs();
	case_2D2t_SGA_fast::glob.gBDelta = f->fgBDelta();
	case_2D2t_SGA_fast::glob.gm = f->fgm();
}

// ***********************************************************************************  //
// *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY ***
// ***********************************************************************************  //


int function(const gsl_vector * x, void *params, gsl_vector * f) {

	struct parameters p;
	p.n = ((Par *) params)->n();
	p.J = ((Par *) params)->J();
	p.t = ((Par *) params)->t();
	p.tPrim = ((Par *) params)->tPrim();
	p.U = ((Par *) params)->U();
	p.beta = ((Par *) params)->beta();
	p.L = ((Par *) params)->L();

	variables v;
	v.barchiAB = gsl_vector_get (x, 0);
	v.chiS = gsl_vector_get (x, 1);
	v.chiT = gsl_vector_get (x, 2);
	v.DeltaS = gsl_vector_get (x, 3);
	v.DeltaT = gsl_vector_get (x, 4);
	v.mu = gsl_vector_get (x, 5);
	v.mAF = gsl_vector_get (x, 6);
	v.d = gsl_vector_get (x, 7);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5),
				gsl_vector_get (x, 6), gsl_vector_get (x, 7));
	}

    // It's start of the calculations

    Equa E(p, v);
    variables w;
    E.run(w);

    // It's end of the calculations

	gsl_vector_set (f, 0, w.barchiAB);
	gsl_vector_set (f, 1, w.chiS);
	gsl_vector_set (f, 2, w.chiT);
	gsl_vector_set (f, 3, w.DeltaS);
	gsl_vector_set (f, 4, w.DeltaT);
	gsl_vector_set (f, 5, w.mu);
	gsl_vector_set (f, 6, w.mAF);
	gsl_vector_set (f, 7, w.d);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5),
				gsl_vector_get (x, 6), gsl_vector_get (x, 7));
		printf("###  %24.15e# %25.15e# %24.15e# %24.15e# %24.15e# %24.15e# %24.15e# %25.15e | E = %24.15e\n",
			gsl_vector_get (f, 0), gsl_vector_get (f, 1), gsl_vector_get (f, 2),
			gsl_vector_get (f, 3), gsl_vector_get (f, 4), gsl_vector_get (f, 5),
			gsl_vector_get (f, 6), gsl_vector_get (f, 7), case_2D2t_SGA_fast::glob.Energy);
		printf("\n");
	}

    return GSL_SUCCESS;
}

// ************* CLASSES ************* //

// Par

void Par::print(FILE *file) const {
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.15e", tabD[i]); }
	for(int i=0; i<NtabI; i++) { fprintf(file, " %6d", tabI[i]); }
}

void Par::print_verboten(FILE *file) const {
	fprintf(file, "PARAMETERS: n=%f, J=%f, t=%f, t'=%f, U=%f, beta=%f, L=%d\n", n(), J(), t(), tPrim(), U(), beta(), L());
}

void Par::print_label(FILE *file) const {
	fprintf(file, "n   J   t   t'   U   beta   L\n");
}

// Var

void Var::print(FILE *file) const{
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.15e", tabD[i]); }
}

void Var::print_verboten(FILE *file) const{
	fprintf(file, "VARIABLES:");
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.15e", tabD[i]); }
	fprintf(file, "\n");
}

void Var::print_label(FILE *file) const {
	fprintf(file, "        barchiAB             chiS                 chiT                    DeltaS                  DeltaT               mu                     mAF                    d");
}

// Res

void Res::set(gsl_multiroot_fsolver *s) {
	for(int i=0; i<n; i++) {
		(*V)[i] = gsl_vector_get (s->x, i);
	}
	this->glob = case_2D2t_SGA_fast::glob;
}

void Res::print(FILE *file) const { V->print(file); fprintf(file, "%21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f\n",
		glob.Energy, glob.gt, glob.gtPrimMax, glob.gtPrimMin, glob.gs, glob.gBDelta, glob.gm, glob.lbarchiAB, glob.lchiS, glob.lchiT, glob.lDeltaS, glob.lDeltaT, glob.ln, glob.lmAF); }
void Res::print_verboten(FILE *file) const { fprintf(file, "Result: "); V->print_verboten(file); fprintf(file, "%21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f\n",
		glob.Energy, glob.gt, glob.gtPrimMax, glob.gtPrimMin, glob.gs, glob.gBDelta, glob.gm, glob.lbarchiAB, glob.lchiS, glob.lchiT, glob.lDeltaS, glob.lDeltaT, glob.ln, glob.lmAF); }
void Res::print_label(FILE *file) const { V->print_label(file); fprintf(file, "                      E                      gt                    gtPrimMax             gtPrimMin            gs                  gBDelta               gm                       lbarchiAB             lchiS                 lchiT             lDeltaS              lDeltaT                  ln                  lmAF\n");}

void example() {

	printf("\nEXAMPLE 2D2t_fast\n");

	Par par(N_PAR);
	par.set_n(0.999);
	par.set_J(1./3.);
	par.set_t(-1.);
	par.set_tPrim(0.);
	par.set_U(12);
	par.set_beta(1500);
	par.set_L(256);
	printf("#"); par.print_verboten(); fflush(stdout);

    //double ZERO = 1e-12;

    Var initV_Z2_U12_n999_B1500(N_VAR);
    initV_Z2_U12_n999_B1500.set_barchiAB(1.693217473268961e-01);
    initV_Z2_U12_n999_B1500.set_chiS(1.442308124273939e-04);
    initV_Z2_U12_n999_B1500.set_chiT(-3.865254347226006e-04);
    initV_Z2_U12_n999_B1500.set_DeltaS(1.669293228458277e-01);
    initV_Z2_U12_n999_B1500.set_DeltaT(3.915289953743466e-05);
    initV_Z2_U12_n999_B1500.set_mu(3.371567247648004e+00);
    initV_Z2_U12_n999_B1500.set_mAF(1.317210700617173e-01);
    initV_Z2_U12_n999_B1500.set_d(1.496803286385485e-02);

    printf("Starting point: "); initV_Z2_U12_n999_B1500.print(); printf("\n\n");
    Res res(N_VAR); res.print_label();
    RootGSL root(par, function, res);
    if(root.Solve(initV_Z2_U12_n999_B1500) == GSL_SUCCESS) { root.print(); }

    printf("Koniec\n");
}

void inf() {
    Par par(N_PAR);
    par.set_n(0.99);
    par.set_J(1./3.);
    par.set_t(-1.);
    par.set_tPrim(0.25);
    par.set_U(8);
    par.set_beta(500.);
    par.set_L(512);
    par.print_verboten();

	// I'm choosing start point
	Var initV(N_VAR);

	// INFINITE
	Res res(N_VAR);
	RootGSL root(par, function, res);
	printf("sdelta      "); res.print_label();
	while(1) {
	    initV.set_barchiAB(2*frand()-1.);
	    initV.set_chiS(2*frand()-1.);
	    initV.set_chiT(2*frand()-1.);
	    initV.set_DeltaS(2*frand()-1.);
	    initV.set_DeltaT(2*frand()-1.);
	    initV.set_mu(16.*frand()-8.);
	    initV.set_mAF(2*frand()-1.);
	    initV.set_d(frand());

		if(root.Solve(initV) == GSL_SUCCESS) {
			printf("       "); root.print(); fflush(stdout);
		}
	}
}

void plot1D() {
	printf("#PLOT1D - 2D2t_fast\n");

	Par initPar(N_PAR);
	initPar.set_n(0.996);
	initPar.set_J(1./3.);
	initPar.set_t(-1.);
	initPar.set_tPrim(0.25);
	initPar.set_U(12);
	initPar.set_beta(1500);
	initPar.set_L(256);

    Var initV_S2(N_VAR);
    initV_S2.set_barchiAB(1.733770448794482e-01);
    initV_S2.set_chiS(6.390883282464091e-04);
    initV_S2.set_chiT(7.001707899582296e-04);
    initV_S2.set_DeltaS(1.643909148624374e-01);
    initV_S2.set_DeltaT(7.833765657089206e-05);
    initV_S2.set_mu(3.489514461003490e+00);
    initV_S2.set_mAF(6.405636335579069e-02);
    initV_S2.set_d(3.134804099861067e-02);

	// I'm choosing start point
	Var initV(N_VAR); initV = initV_S2;

	printf("#dochodze do punktu startowego\n");
    ////automatic_run_U(initPar, initV, 10., 9.98, 0.01); printf("\n");
	//automatic_run_n(initPar, initV, 0.996, 0.9998, 0.0001); printf("\n");
    //automatic_run_tPrim(initPar, initV, 0.25, 0.20, 0.01); printf("\n");
    //automatic_run_U(initPar, initV, 12., 10.4, 0.2); printf("\n");
    //automatic_run_J(initPar, initV, 0.33, 0.10, 0.01); printf("\n");
	//automatic_run_beta(initPar, initV, 1500, 500, 100); printf("\n");



	printf("#Jestem w punkcie startowym\n");

	Par oldPar = initPar;
    Var oldVar = initV;

    automatic_run_tPrim(initPar, initV, 0.25, 0, 0.005); printf("(chyba niepotrzebne)\n"); initPar=oldPar; initV=oldVar;

}

void plot2D() {
	printf("#PLOT2D - 2D2t_fast\n");

	Par initPar(N_PAR);
	initPar.set_n(0.996);
	initPar.set_J(1./3.);
	initPar.set_t(-1.);
	initPar.set_tPrim(0.25);
	initPar.set_U(12);
	initPar.set_beta(1500);
	initPar.set_L(0);


    // STANDARD - t'=0.25, t''=0
    //PARAMETERS: n=0.996000, J=0.333333, t=-1.000000, t'=0.250000, t''=0.000000, U=12.000000, beta=500.000000, L=512
    //                      barchiAB             simchiAB             chiS                 chiT                     chiSPrim             chiTPrim            DeltaS                  DeltaT                  simDeltaB               mu                     mAF                    d                      E                      gt                    gtPrimMax             gtPrimMin             gtPrim2Max            gtPrim2Min            gs                  gBDelta               gm                 lbarchiAB                 lsimchiAB             lchiS                 lchiT                lchiSPrim              lchiTPrim             lDeltaS              lDeltaT              lsimDeltaB            ln                  lmAF
    // DS         1.736126774702611e-01 -4.374535151143151e-23 -7.105665884589442e-04  3.253559670497917e-13  2.370085064596791e-03  4.210151247296573e-12  1.647364994545274e-01  2.435357064581637e-13  3.054348850582745e-21 3.504606218225224e+00 -3.102341566783024e-11  3.143224122455469e-02   -0.243329264265046     0.020707478060625     0.020707478060137     0.020707478061114     0.020707478060137     0.020707478061114     3.952461554088670     0.020707478060626     1.988079866124264     0.192256836311526    -0.000000000000000    -0.005176869515156     0.000000000000000    -0.000000000000000     0.000000000000000     0.162778670162292    -0.000000000000080    -0.000000000000000    -7.027222728294262    -0.000000000025031
//S2// DS+DT+m    1.733770448794482e-01  8.150135383561053e-20 -6.390883282464091e-04  7.001707899582296e-04  2.172872463005137e-03  7.516737988623041e-03  1.643909148624374e-01 -7.833765657089206e-05  1.379873610715300e-19 3.489514461003490e+00 -6.405636335579069e-02  3.134804099861067e-02   -0.243338270563175     0.020646155998253     0.019662556713935     0.021678958830522     0.019662556713935     0.021678958830522     3.920307359911151     0.020670757772228     1.979976605899966     0.190568982268390     0.000000000000000    -0.005167689443057     0.000000000000000    -0.000000000000000     0.000000000000000     0.161115728359435     0.000025592307634    -0.000000000000000    -6.997883807312957    -0.050805948849911

    Var initV_S2(N_VAR);
    initV_S2.set_barchiAB(1.733770448794482e-01);
    initV_S2.set_chiS(6.390883282464091e-04);
    initV_S2.set_chiT(7.001707899582296e-04);
    initV_S2.set_DeltaS(1.643909148624374e-01);
    initV_S2.set_DeltaT(7.833765657089206e-05);
    initV_S2.set_mu(3.489514461003490e+00);
    initV_S2.set_mAF(6.405636335579069e-02);
    initV_S2.set_d(3.134804099861067e-02);

	/*// U20, n09999, tPrim025, J13, B17
    //17.000000  1.673883323527769e-01 -7.216007878910981e-06 -2.247544589671341e-05 1.672485489116952e-01 2.464253195574323e-06
    //2.670567575333385e+00 1.248903527436433e-01 1.825637177834200e-03   -0.230920603594082     0.000285433655713     0.000339727829034     0.000239816596848     3.877251870157079     0.000289772212941     1.969073861021236     0.162537114820033    -0.000072443053235     0.000000000000000     0.162116187262232    -0.000000796210859    -5.341432804489117     0.099248648927394    Var initV_S2_B17(N_VAR);
	initPar.set_n(0.9999);
	initPar.set_J(1./3.);
	initPar.set_t(-1.);
	initPar.set_tPrim(0.25);
	initPar.set_U(20);
	initPar.set_beta(17);
	initPar.set_L(1024);
    Var initV_S2_B17(N_VAR);
    initV_S2_B17.set_barchiAB(1.673883323527769e-01);
    initV_S2_B17.set_chiS(-7.216007878910981e-06);
    initV_S2_B17.set_chiT(-2.247544589671341e-05);
    initV_S2_B17.set_DeltaS(1.672485489116952e-01);
    initV_S2_B17.set_DeltaT(2.464253195574323e-06);
    initV_S2_B17.set_mu(2.670567575333385e+00);
    initV_S2_B17.set_mAF(1.248903527436433e-01);
    initV_S2_B17.set_d(1.825637177834200e-03);
*/

	// I'm choosing start point
	Var initV(N_VAR); initV = initV_S2;

	initPar.set_L(760);
	printf("#dochodze do punktu startowego\n");
    ////automatic_run_U(initPar, initV, 10., 9.98, 0.01); printf("\n");
	automatic_run_n(initPar, initV, 0.996, 0.9998, 0.0001); printf("\n");
		automatic_run_n(initPar, initV, 0.9998, 0.99998, 0.00001); printf("\n");
		automatic_run_n(initPar, initV, 0.99998, 0.99999, 0.000001); printf("\n");
    automatic_run_tPrim(initPar, initV, 0.25, 0.20, 0.01); printf("\n");
    automatic_run_U(initPar, initV, 12., 10.4, 0.2); printf("\n");
    automatic_run_U(initPar, initV, 10.4, 10.0, 0.01); printf("\n");
    //automatic_run_J(initPar, initV, 0.33, 0.10, 0.01); printf("\n");
    	//automatic_run_J(initPar, initV, 0.10, 0.01, 0.001); printf("\n");

		//automatic_run_beta(initPar, initV, 1500, 500, 100); printf("\n");
		//automatic_run_beta(initPar, initV, 500, 200, 50); printf("\n");
		//automatic_run_beta(initPar, initV, 200, 100, 10); printf("\n");
		//automatic_run_beta(initPar, initV, 100, 50, 5); printf("\n");
		//automatic_run_beta(initPar, initV, 50, 20, 1); printf("\n");
		//automatic_run_beta(initPar, initV, 20, 15, 0.2); printf("\n");


	/*
	 automatic_run_beta(initPar, initV, 17, 10, 1); printf("\n");
	*/
	printf("#Jestem w punkcie startowym\n");

	Par oldPar = initPar;
    Var oldVar = initV;

    //printf("#Zaczynam 4-ta cwiartke\n");	automatic_2Drun_Un(initPar, initV, 20., 9., 0.5, 0.9999, 0.994, 0.0001);		printf("#Czwarta cwiartka skonczona\n"); initPar=oldPar; initV=oldVar;
    //printf("#Zaczynam 2-ga cwiartke\n");	automatic_2Drun_nU(initPar, initV, 0.9999, 0.994, 0.0003, 15., 9., 0.1);		printf("#Druga cwiartka skonczona\n"); initPar=oldPar; initV=oldVar;

    //printf("#Zaczynam 1-szą cwiartke\n");	automatic_2Drun_nU(initPar, initV, 0.999, 1., 0.00005, 12., 5., 0.1);		printf("#Pierwsza cwiartka skonczona\n"); initPar=oldPar; initV=oldVar;
    //printf("#Zaczynam 2-ga cwiartke\n");	automatic_2Drun_nU(initPar, initV, 0.999, 0.998, 0.00005, 12., 5., 0.1);	printf("#Druga cwiartka skonczona\n"); initPar=oldPar; initV=oldVar;
    //printf("#Zaczynam 3-cia cwiartke\n");	automatic_2Drun_nU(initPar, initV, 0.999, 1., 0.00005, 12., 30., 0.5);		printf("#Trzecia cwiartka skonczona\n"); initPar=oldPar; initV=oldVar;
	//printf("#Zaczynam 4-ta cwiartke\n");	automatic_2Drun_nU(initPar, initV, 0.999, 0.99, 0.0001, 12., 30., 0.5);		printf("#Czwarta cwiartka skonczona\n"); initPar=oldPar; initV=oldVar;

    //printf("#W dol\n");  automatic_2Drun_JU(initPar, initV, 0.01, 0., 0.001, 14., 11., 0.005);		printf("#Koniec (w dol)\n"); initPar=oldPar; initV=oldVar;
    //printf("#W gore\n"); automatic_2Drun_JU(initPar, initV, 0.30, 3., 0.05, 15., 5., 0.1);		printf("#Koniec (w gore)\n"); initPar=oldPar; initV=oldVar;

    //printf("#W dol\n");  automatic_2Drun_Jn(initPar, initV, 0.10, 0., 0.001, 0.99999, 0.998, 0.00001);		printf("#Koniec (w dol)\n"); initPar=oldPar; initV=oldVar;
    //printf("#W gore\n"); automatic_2Drun_Jn(initPar, initV, 0.30, 3., 0.05, 0.999, 0.98, 0.0002);		printf("#Koniec (w gore)\n"); initPar=oldPar; initV=oldVar;


    //automatic_run_n(initPar, initV, 0.998, 0.99, 0.0001); printf("\n");
    //automatic_run_nlog(initPar, initV, 0.9999, 1., 150); printf("\n");

    //automatic_run_n(initPar, initV, 0.999, 1, 0.00005); printf("\n"); initPar=oldPar; initV=oldVar;
    //automatic_run_n(initPar, initV, 0.999, 0.99, 0.0001); printf("\n");
    //automatic_run_n(initPar, initV, 0.99, 0.5, 0.002); printf("\n");

    //printf("#W dol\n");  automatic_2Drun_tPrimn(initPar, initV, 0.5, 0., 0.05, 0.995, 0.9941, 0.000001);		printf("#Koniec (w dol)\n"); initPar=oldPar; initV=oldVar;
    //automatic_run_n(initPar, initV, 0.995, 0.9954, 0.0001); printf("\n");
    //printf("#W gore\n"); automatic_2Drun_tPrimn(initPar, initV, 0.5, 1., 0.05, 0.9954, 0.9946, 0.00001);		printf("#Koniec (w gore)\n"); initPar=oldPar; initV=oldVar;

    //printf("#W dol\n");  automatic_2Drun_tPrimU(initPar, initV, 0.25, 0., 0.05, 10.0, 9.8, 0.0005);		printf("#Koniec (w dol)\n"); initPar=oldPar; initV=oldVar;
    //printf("#W gore\n"); automatic_2Drun_tPrimU(initPar, initV, 0.25, 1., 0.05, 10.0, 9.8, 0.0005);		printf("#Koniec (w gore)\n"); initPar=oldPar; initV=oldVar;

    automatic_run_U(initPar, initV, 10.0, 9.50, 0.0005); printf("\n");

    printf("#KONIEC\n");

    return;
}




/////////////////////////
///// AUTOMATIC PLOTS ///
/////////////////////////

void automatic_2Drun_tPrimU(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_tPrim(start1);
	par.set_U(start2);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#tPrim    U        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_tPrim(start1 + ((double)i)*step1);
		par.set_U(start2);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_U(start2 + ((double)j)*step2);
			printf("%f %f   ", par.tPrim(), par.U());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_tPrim(stop1); par.set_U(stop2);
	printf("\n#END\n");
	fflush(stdout);

	return;
}

void automatic_2Drun_tPrimn(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_tPrim(start1);
	par.set_n(start2);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#tPrim    n        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_tPrim(start1 + ((double)i)*step1);
		par.set_n(start2);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_n(start2 + ((double)j)*step2);
			printf("%f %f   ", par.tPrim(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_tPrim(stop1); par.set_n(stop2);
	printf("\n#END\n");
	fflush(stdout);

	return;
}

void automatic_2Drun_Jn(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_J(start1);
	par.set_n(start2);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#J        n        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_J(start1 + ((double)i)*step1);
		par.set_n(start2);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_n(start2 + ((double)j)*step2);
			printf("%f %f   ", par.J(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_J(stop1); par.set_n(stop2);
	printf("\n#END\n");
	fflush(stdout);

	return;
}

void automatic_2Drun_JU(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_J(start1);
	par.set_U(start2);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#J        U        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_J(start1 + ((double)i)*step1);
		par.set_U(start2);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_U(start2 + ((double)j)*step2);
			printf("%f %f   ", par.J(), par.U());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_J(stop1); par.set_U(stop2);
	printf("\n#END\n");
	fflush(stdout);

	return;
}

void automatic_2Drun_Un(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_U(start1);
	par.set_n(start2);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#U        sdelta  "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_U(start1 + ((double)i)*step1);
		par.set_n(start2);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_n(start2 + ((double)j)*step2);
			printf("%f %f   ", par.U(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_U(stop1); par.set_n(stop2);
	printf("\n#END\n");
	fflush(stdout);

	return;
}

void automatic_2Drun_nU(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_n(start1);
	par.set_U(start2);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta  U        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_n(start1 + ((double)i)*step1);
		par.set_U(start2);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_U(start2 + ((double)j)*step2);
			printf("%f %f   ", par.n(), par.U());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_n(stop1); par.set_U(stop2);
	printf("\n#END\n");
	fflush(stdout);

	return;
}

void automatic_run_nlog(Par &inpar, Var &invar, double start, double stop, double loop) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_n, start %f, stop %f, loop %f | ", start, stop, loop); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#n                 "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = loop;
	double INTERVAL = stop-start;
	for(int i=0; i<LOOPS; i++) {
		inpar.set_n(start + INTERVAL*tanh((double)i/10.));
		printf("%.15f ", inpar.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	inpar.set_n(stop);
	invar = initV;

	printf("#end automatic_run_n\n\n");
	fflush(stdout);
	return;
}

void automatic_run_n(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_n, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#n           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step+0.5);
	if(LOOPS < 0) { LOOPS = -LOOPS+1; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_n(start + ((double)i)*step);
		printf("%f ", inpar.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	inpar.set_n(stop);
	invar = initV;

	printf("#end automatic_run_n\n\n");
	fflush(stdout);
	return;
}

void automatic_run_tPrim(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_tPrim, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#tPrim       "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step+0.5);
	if(LOOPS < 0) { LOOPS = -LOOPS+1; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_tPrim(start + ((double)i)*step);
		printf("%f ", inpar.tPrim());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_tPrim\n\n");
	fflush(stdout);

	inpar.set_tPrim(stop);
	invar = initV;
	return;
}

void automatic_run_J(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_J, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#J           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step+0.5);
	if(LOOPS < 0) { LOOPS = -LOOPS+1; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_J(start + ((double)i)*step);
		printf("%f ", inpar.J());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_J\n\n");
	fflush(stdout);

	inpar.set_J(stop);
	invar = initV;
	return;
}

void automatic_run_U(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_U, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#U           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step+0.5);
	if(LOOPS < 0) { LOOPS = -LOOPS+1; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_U(start + ((double)i)*step);
		printf("%f ", inpar.U());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_U\n\n");
	fflush(stdout);

	inpar.set_U(stop);
	invar = initV;
	return;
}

void automatic_run_beta(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_beta, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#beta        "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step+0.5);
	if(LOOPS < 0) { LOOPS = -LOOPS+1; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_beta(start + ((double)i)*step);
		printf("%f ", inpar.beta());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_beta\n\n");
	fflush(stdout);

	inpar.set_beta(stop);
	invar = initV;
	return;
}

void automatic_run_L(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_L, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#L           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step+0.5);
	if(LOOPS < 0) { LOOPS = -LOOPS+1; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_L(start + ((double)i)*step);
		printf("%d ", inpar.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_L\n\n");
	fflush(stdout);

	inpar.set_L(stop);
	invar = initV;
	return;
}





} // end of namespace case2D2t_SGA_tPrim0
